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ABSTRACT 



A detailed parametric study of the variables critical to a numerical 
model of the interaction of a free surface with an ascending vortex 
pair has been carried out. The model is based on a discrete vortex 
representation of the free surface and the use of the appropriate 
boundary conditions. A study was performed to minimize the numeri- 
cal Helmholtz instabilities inherent in the numerical model for lower 
Froude numbers. A method of free-surface re-discretization was 
developed and implemented. Additionally, an investigation of the 
motion of a vortex pair approaching and passing through a sharp 
salt/freshwater density interface was conducted. Experiments were 
performed in a large basin and the vortices were generated with two 
counter-rotating plates. The rise and intrusion of the Kelvin oval into 
the upper layer of lower density fluid was recorded on a video tape and 
then analyzed through the use of a Motion Analysis system. 
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I. INTRODUCTION 



The understanding of the temporal development of the surface 
signatures resulting from the interaction of an ascending inclined vor- 
tex pair with the free surface is of special importance. When a trailing 
vortex pair, generated by a lifting surface, approaches the free surface, 
with or without mutual induction instability and/or vortex breakdown, 
the vortices and/or the crude vortex rings give rise to surface distur- 
bances (scars and striations). The striations are essentially three- 
dimensional free-surface disturbances, normal to the direction of the 
motion of the lifting surface, and come into existence when the vortex 
pair is at a distance equal to about one initial vortex separation from 
the free surface. The scars are small free-surface depressions on 
either side of a nearly symmetrical hump, directly above the plane of 
symmetry of the vortex pair. They come into existence towards the 
end of the pure striation phase and when the vortices are at a distance 
equal to about 60 percent of the initial vortex separation from the free 
surface. The scars are created and driven by the vortices and both 
move outward. Subsequently, the vortices diffuse rapidly and their 
circulation decreases slowly at first (prior to the occurrence of scars) 
and rapidly thereafter due to the counter-sign vorticity generated in 
the overlapping regions of the vortex pair in the recirculation cell and 
at the free surface. 
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The surface disturbances or the local roughening of the free sur- 
face, resulting from the vortex motion produced by submerged foils or 
counter-rotating propellers, are visible because the particle motions 
they create as they propagate produce a straining field at the surface 
which modifies the physical characteristics of the sea surface 
(reflection, refraction, and the scattering of sound, heat, and light). In 
fact, in many photographs of the sea surface taken from aircraft or 
earth resources technology satellites, using passive or active sensors, 
the differential roughening of the free surface can be seen clearly. 

Previous studies of a vortex (or vortex pair) in ground- or free- 
surface effect have been generally concerned with the boundary layer 
development, generation of oppositely-signed vorticity, the flow 
resulting from a single vortex held fixed in a uniform flow and the so- 
called “rebounding” phenomenon, according to which the vortices 
move at first toward the boundary, as on an inviscid trajectory, and 
then away from it. The boundary surface may be either a rigid bound- 
ary at which the no-slip condition is satisfied or a free surface corre- 
sponding to zero shear stress. In either case, the boundary surface 
(sometimes erroneously referred to as a vorticity sink) gives rise to 
oppositely signed vorticity and the net circulation of the vortices 
decreases with time. 

The interaction of vortices with a fluid interface has been under 
investigation at the Naval Postgraduate School by Sarpkaya and his stu- 
dents for the past several years [Refs. 1-11]. The emphasis of these 
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investigations has been the identification of the specific mechanisms 
which lead to the demise of the vortices and the prediction, through 
experiments and analysis, of the motion of vortices and the form of the 
free surface. 

Sarpkaya and Henderson’s [Ref. 9] theoretical analysis of the scars 
created by the trailing vortices was carried out assuming the vortices 
to be two-dimensional and the free surface to be a rigid plane. For 
small Froude numbers (Fr = V Q /^J gb 0 ) , (where V Q is the initial 
mutual induction velocity of the vortices and b 0 is the initial vortex 
separation), each scar front was shown to coincide with the stagnation 
point on the Kelvin oval, formed by one of the pair of the trailing 
vortices and its image. For Froude numbers larger than about 0.15, 
not only the deformation of the free surface but also the nonlinear 
interaction between the said deformation and the motion of the 
vortices is significant. Thus, the free surface may no longer be 
assumed rigid. However, the analysis of a three-dimensional 
interaction between the surface and the inclined vortex pair is 
extremely difficult. It is in view of this fact that it was decided to 
generate nearly two-dimensional vortex pairs at relatively high Froude 
numbers and to carry out a numerical analysis assuming a deformable 
free surface. Most recently, Leeker [Ref. 11] continued the work on 
the interaction of a vortex pair with the free surface. His task was to 
improve upon the discrete vortex model, originally constructed by 
Elnitsky [Ref 10], by reducing or eliminating the numerical 
instabilities present at lower Froude numbers. This improvement 



3 



would enable the model to more closely represent the results 
observed in the laboratory by allowing the calculations to progress 
further in time. 



The present investigation has a two-fold purpose: 

1. To conduct a detailed parametric study of the variables affecting 
the performance of the discrete vortex model and identify the 
variables (if any) that adversely or positively affect the perfor- 
mance of the model with particular emphasis on their effects on 
observed numerical instabilities. Once it has been determined 
how sensitive the numerical model is to a variety of variables, 
efforts are made to attempt a qualitative optimization of these 
variables under specific conditions. Armed with the information 
gained in the parametric study, additional efforts are made to 
obtain methods that remove or neutralize numerical instabilities 
present at lower Froude numbers. Several smoothing techniques 
and re-discretization schemes are attempted. 

2. To investigate the motion of a vortex pair approaching and 
intruding through a sharp density interface. Utilizing an 
improved test apparatus and a sharp salt/fresh water interface, a 
complete investigation was conducted to analyze the interaction 
of the vortices and the Kelvin oval with the interface. The motion 
of the vortex centers was recorded on video tape and analyzed 
through the use of a motion analysis system. 
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H. PARAMETRIC STUDY OF THE NUMERICAL MODEL 



A INTRODUCTION 

The discrete vortex model used in this investigation and in the 
previous work is discussed by Elnitsky and Leeker in detail [Refs. 10, 
11]. In summary, discrete vortices are placed along the free surface to 
form a vortex sheet. The inter- vortex spacing may vary as desired, 
emanating from a pre-selected center of vortex distribution along the 
free surface. After the program is executed, each individual surface 
vortex is convected utilizing the velocities imparted to it by all other 
vortices and the main rising vortex pair. Execution continues until 
numerical Helmholtz instabilities develop and chaos ensues. 

The complex function representing all of the vortices in the sys- 
tem is given by: 



The strengths of the free-surface vortices are normalized by the 
strength of the main vortex. Spatial units are normalized by the dis- 
tance b 0 , the initial spacing between the two main vortex centers. 
Velocities are normalized by V 0 =r o /2b 0 and time is normalized by 



The normalized boundary condition at the free surface is then 
given by: 



ln(z -z 0 ) _ ln(z + z D ) 
2i n 2in 



i r>1 X~ Zk) + i rtln ^; z>) m 

k =i 2in k=i 2lK 
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and the induced velocity becomes: 



dw 

dt 



u - IV = — 



+ 



+ 



ir. 



-I- 



ir \ 



2n(z-z 0 ) 2n(z+z 0 ) k=1 2 7t(z-z t ) ksl 2]t(z + z ( ) 



(3) 



There are other numerical methods available to solve the problem 
of the interaction of a free surface with a vortex pair, and Elnitsky dis- 
cusses them in great detail [Ref 10]. The principal advantages of the 
Boundary Integral Equation method utilized in this and preceding 
investigations are the ability to discretize the free surface without a 
loss of accuracy and the ability to solve the boundary conditions 
exactly. Another advantage that this method enjoys is the relative 
computer efficiency with which it performs with respect to total com- 
putational time. However, in order to utilize this model to accurately 
predict the behavior of the motion of a vortex pair in nature, it is nec- 
essary that a means be discovered to minimize the impact of the 
numerical instabilities. Finite Element and Finite Difference tech- 
niques were also considered as a means of solving this problem but 
xvere determined to be too cumbersome and computationally 
intensive. 

The initial step in the task of minimizing numerical instabilities 
within this distributed vortex model is to identify which variables most 
affect the performance of the numerical code. In concert with this 
task is the requirement to qualitatively assess the relative sensitivity of 



the code to each changing variable. Finally, after each variable has 
been qualitatively optimized, additional methods must be imple- 
mented to further improve the numerical stability of the model for 
lower Froude numbers. The effect of the following parameters were 
investigated: 

5: desingularization parameter 

N: total number of free surface discrete vortices 

xmax: maximum x/b 0 position 

xv: center of the surface vortex distribution (x/b 0 ) 

Upon the completion of the parametric study, a re-discretization 
method was devised and implemented. 

R DESINGULARIZATION PARAMETER (5) 

1. Introduction 

Leeker [Ref 11] discussed a desingularization method intro- 
duced by Rosenhead [Ref 12] which uses a smoothing parameter to 
remove the singularity at the center of a given vortex and to minimize 
the induced velocity of a vortex on its nearest neighbors. Utilizing this 
method of desingularization, equation (3) becomes: 



o 
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Clearly, this method of desingularization eliminates the 
unnatural infinite velocity at the vortex center and allows the operator 
to more closely space individual vortices at the free surface. Reduced 
inter-vortex spacing minimizes the leakage problem arising from 
widely spaced vortices, yet without the desingularization parameter, 
this close spacing would immediately lead to individual surface vor- 
tices rotating about one another, creating irreparable numerical insta- 
bilities. It is clear that as the spacing between vortex elements is 
increased, the effect of 8 on the performance of the code would 
diminish accordingly. In order to more fully understand the impact of 
this de-singularization, a systematic parametric study was undertaken. 

From earlier investigations by Elnitsky and Leeker [Refs. 10, 
11] the optimum parameters for executing the numerical model and 
ensuring that it would progress an adequate amount of time were 
qualitatively determined to be: 



Number of surface vortices (N) 


76 


Froude number (Fr) 


.35 


Starting depth (d 0 ) 


-3.0 


Starting time (T*) 


-3.00 


2 . Application 





The smoothing parameter 8 was varied from 0.05 to 0.15 and 
the program w r as executed from a non-dimensional starting time of 
-3.00 The choice of 8 was governed by the number of surface ele- 
ments that would be affected at one time by a given vortex. If 8 is cho- 
sen too large with respect to the quantity (z-Zi), the program will 
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essentially ignore several consecutive vortex contributions in its cal- 
culation of the velocity field; this is particularly prevalent at the early 
stages of the motion where the vortex spacing is comparatively small 
near the axis of symmetry. Ideally, the goal was to avoid having 5 mask 
more than two consecutive surface vortices at the beginning of the 
distribution and yet purposely mask vortices where the numerical 
instabilities form. This allows the calculation to proceed with each 
vortex having only a minor effect on its nearest neighbors but still 
having the full effect on more distant vortices. As a result, the opera- 
tor may space vortices very close together and yet still minimize the 
instabilities which would otherwise occur. 

Table 1 summarizes a representative number of executions of 
the numerical model while holding all other variables constant and 
varying 5 from .05 to .15. The breakdown time represents the last 
reliable time (± five time steps) where the free surface is not unnatu- 
rally formed. As can be seen clearly from Figure 1, the lowest value of 
8 (.05) is unacceptable. The free surface deforms in an entirely unnat- 
ural fashion at all non-dimensional times later than -1.4. Run 2 
demonstrates significant improvement with only a very minor increase 
in 8 (Figure 2). For subsequent runs with a slightly increasing 8 value, 
the gain in extended execution time is minimized and appears to con- 
verge somewhat to a non-dimensional time of about -1.0. Figures 3 
through 7 illustrate the form of the free surface for each run at the last 
reliable time just prior to the halt of the execution due to numerical 
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Figure 1. Run 1, Free Surface Deformation, 5 = .05, T* = -1.4 
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Figure 2. Run 2, Free Surface Deformation, 8 = .07, T* = -1.25 




Figure 3. Run 3, Free Surface Deformation, 5 = .09, T* = -1.1 
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Figure 4. Run 4, Free Surface Deformation, 5 = 0.1, T* = -1.1 
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Figure 5. Run 5, Free Surface Deformation, 5 = .11, T* = -1.05 
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Figure 6 
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. Run 6, Free Surface Deformation, 5 = .13, T* = -1.05 
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Figure 7. Run 7, Free Surface Deformation, 5 = .15, T* = -1.0 
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instabilities. Essentially, any increase in 6 above approximately 0. 1 
serves no other purpose than to mask a number of surface vortices 
unnecessarily, with only a trivial gain in additional time steps of 
execution. Figures 8 through 14 represent the free-surface deforma- 
tion at execution-halt time and clearly demonstrate that the gain in 
program execution time is minimal. The numerical model is surpris- 
ingly insensitive to this parameter once it has accomplished its pur- 
pose of removing the singularity. The qualitatively chosen optimum 
value of 5 was determined to be 0.1. 

TABLE 1 

VARIATION OF 5 VERSUS BREAKDOWN TIME 

Breakdown 



Run# 


5 


Time (T 


1 


.05 


-1.4 


2 


.07 


-1.25 


3 


.09 


-1.1 


4 


.1 


-1.1 


5 


.11 


-1.05 


6 


.13 


-1.05 


7 


.15 


-1.0 


NUMBER OF FREE SURFACE VORTICES (N) 



1. Introduction 

The number of free surface vortices present during any par- 
ticular run is automatically calculated by a subroutine within the 
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Figure 8. Run 1, Free Surface Deformation, 5 = .05, 
Execution Halt T*=-1.20 
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Figure 9. Run 2, Free Surface Deformation, 8 = .07, 
Execution Halt T* = -1.10 
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Figure 10. Run 3, Free Surface Deformation, 6 = .09, 
Execution Halt T* = -1.05 
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Figure 11. Run 4, Free Surface Deformation, 8 = .10, 
Execution Halt T* = -1.00 
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Figure 12. Run 5, Free Surface Deformation, 8 = .11, 
Execution Halt T* = -1.00 
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Figure 13. Run 6, Free Surface Deformation, 5 = .13, 
Execution Halt T* = -0.95 
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Figure 14. Run 7, Free Surface Deformation, 8 = .15, 
Execution Halt T* = -0.95 
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program when given certain desired parameters. The operator 
directs the program to place the surface vortices after deciding on the 
following conditions: 

a The initial distance between the first free surface vortex and the 
centerline. 

b. The position (x/bo) where the geometric progressional place- 
ment of vortices ends and constant spacing begins. 

c. The actual distance (x/bo) between vortices after the point in (b) 
above has been reached. 

The details of the calculation of vortex spacing is not particu- 
larly important but the actual decision of how many vortices to place 
along the free surface will critically affect the performance of the 
model. As discussed earlier, as vortex density increases, although 
leakage decreases (a desirable result), the velocity impact of one vor- 
tex upon another increases and rapidly forces the breakdown of the 
model prematurely. 

2. Application 

Numerous program executions were completed to determine 
the sensitivity of the numerical model to the number of discrete free 
surface vortices. Table 2 represents the results of a portion of these 
executions. 

The desingularization parameter (5) was adjusted in order to 
maintain the criteria established earlier (avoid masking vortices 
unnecessarily yet still “cover" the breakdown point). This adjustment 
of 5 was necessary to ensure that the model remained as close to its 
natural form as possible. For example, if 6 was maintained at 0.1 
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during run #10, it would effectively mask a much larger number of 
discrete vortices and thus introduce an unnecessary artificiality into 
the execution of the model. Similarly, if 6 w r as maintained at 0.1 
during run #8, the singularity in the velocity expression would barely 
even be removed. Figures 15 through 17 depict the formation of the 
free surface immediately before breakdown. 

TABLE 2 

VARIATION OF N VERSUS BREAKDOWN TIME 

Breakdown 



Run 


5 


N 


Time (T*) 


8 


0.18 


38 


-1.0 


9 


0.10 


76 


-1.1 


10 


0.085 


149 


-1.1 



The results of these executions reveal the model’s insensitiv- 
ity to a variable number of discrete free surface vortices. The only 
concern remaining to the operator as a result of this particular study is 
to balance the leakage versus the propensity of the vortices to orbit 
about each other, as discussed earlier, and to consider the total com- 
putational time required to complete each run. For example. Runs #8 
and #9 require only a moderate amount of time to complete, but Run 
#10 required over 90 minutes of real-time computation. As a result, 
any execution with N > 149 must be executed during low-computer- 
use time only. Qualitatively weighing these factors, it was determined 
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Figure 15. Run 8 , Free Surface Deformation, N 
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Figure 16. Run 9, Free Surface Deformation, N = 76 
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Figure 17. Run 10, Free Surface Deformation, N = 149 
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that the optimum number of free surface vortices was 76. Average 
vortex density (N/xmax) was held constant throughout the remainder 
of the parametric study. 

D. MAXIMUM X/B 0 VALUE (XMAX) 

1. Introduction 

Very closely related to the total number of vortices in the 
model is the total length of free surface that they will be responsible 
for representing. It is not desirable to have the vortices represent an 
extensive stretch of free surface since it would prolong the total com- 
putation time and serve no meaningful function. It is also undesirable 
to restrict the extent of the free surface so severely that it would not 
truly represent nature as closely as possible. The simple goal of the 
study of the numerical model’s sensitivity to a changing xmax value is 
to confirm the theory that the model would behave similarly at differ- 
ent N, provided the vortex density remained relatively constant. 

2. Application 

Two executions were completed in an attempt to confirm the 
theory of vortex density within the model. Table 3 represents the 
results of these executions. 
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TABLE 3 



VARIATION OF XMAX VERSUS BREAKDOWN TIME 

Breakdown 

Run # 5 N xmax Time (T*) 

11 .1 76 10 -1.0 

12 .1 38 5 -1.1 

From Figures 18 and 19, the form of the free surface demon- 
strated no significant change when xmax was decreased from 10 x/b 0 

to 5 x/b 0 . Essentially, there is no major advantage in using a value 
exceeding about 5 x/b 0 since all vortex action occurs well within that 

value. The only remaining advantage of varying this parameter is to 
increase or decrease real CPU time. Obviously, as the value of xmax 

increases, computational time increases as well. For the remainder of 
this parametric study, a value of 5 x/b 0 was utilized. 

E. CENTER OF THE SURFACE VORTEX DISTRIBUTION (XV) 

1. Introduction 

Elnitsky (Ref. 10] pointed out that a major advantage of the 
Boundary Integral method is the propensity of the free surface vortices 
to migrate to the regions where the radius of curvature of the free 
surface is smallest. The surface vortices “roll down” the hump of the 
deformed free surface and collect in the region of the scar. Thus one 
achieves the greatest discretization at the positions along the free 
surface where it is most needed. 
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Figure 18. Run 11, Free Surface Deformation, xmaz * 10 
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Figure 19. Run 12, Free Surface Deformation, mar = 5 
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As discussed earlier, the center of free surface vortex distri- 
bution is an operator-selected parameter. The algorithm that calcu- 
lates the initial position of each free surface vortex will allow the 
distribution of vortices to emanate from any x/b 0 position. From that 
selected position, a subroutine will geometrically place vortices to the 
left until the centerline (x/b 0 = 0) has been reached and to the right 
until the x/b 0 position of xmax has been reached. Since the vortex 
density along the free surface is greatest at the center of distribution, 
by moving this point the operator can affect the vortex density at any 
point along the free surface later in time. Naturally, as one progresses 
away from the center of distribution, this effect is minimized and is 
eventually eliminated completely. The purpose of the variation of this 
parameter is to ascertain whether there is an ideal (xv) position 
enabling the surface vortices to “roll” into an optimum position at the 
correct time and minimize or delay the onset of the numerical insta- 
bilities that are introduced into the free surface. 

2. Application 

Runs 13 through 15 represent a sample of the executions of 
the code completed in order to identify the sensitivity of the model to 
the parameter (xv). Table 4 lists the details of each run and its effect 
on completion time. 
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TABLE 4 



VARIATION OF xv VERSUS BREAKDOWN TIME 


Run# 


xv (x/bo) 


Breakdown 
Time (T*) 


13 


1.0 


-1.10 


14 


0.8 


-1.05 


15 


0.5 


-1.25 



Although the breakdown time of runs 13 through 15 did not 
vary significantly from earlier runs, the variation of xv resulted in 
totally unacceptable free surface formations. From Figures 20 through 
22, it is noticed that the flat shape of the hump later in time is caused 
by the relatively low vortex density Also noted is the fact that the free 
surface vortices continue to “roll downhill” but terminate farther than 
in previous executions. In fact, this tendency of the surface vortices to 
roll farther away from centerline actually contributes to the formation 
of the numerical instabilities by creating sharp corners and large 
spaces between individual vortices. Accordingly, it was decided to 
maintain the center of vortex distribution at the centerline. 
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Figure 20. Run 13, Free Surface Deformation, xv = 1.0 
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Figure 21. Run 14, Free Surface Deformation, xv = 0.8 
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Figure 22. Run 15, Free Surface Deformation, xv = 0.5 
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III. SMOOTHING 



A INTRODUCTION 

Despite all efforts to minimize the effect, or delay the onset, of 
the numerical Helmholtz instabilities during the parametric study, 
they consistently occur and are clearly inevitable. Extensive efforts 
were expended to discover an additional parameter within the model 
that would minimize the effect of the instabilities or delay their 
formation. During the course of the parametric study, as a path of 

investigation, it was decided to observe the variation of the circulation 
strength of the free surface vortices with respect to increasing x/b 0 

distance as measured from center line. It was theorized that sharp 
differences in circulation values between adjacent vortices were 
responsible for the tendency of the adjacent vortices to rotate about 
each another. 

As a representative sample, run 12 was utilized to depict the vari- 
ation of the circulation with respect to x/b 0 . Figure 23 clearly demon- 
strates the so called “sawtooth” instabilities. Each point on the figure 
represents an individual free surface vortex circulation strength. 

It is useful to observe the interaction of the Helmholtz instabilities 
with the variation of the circulation values. Figure 24 depicts the free 
surface vortices that eventually (within one or two time steps) will lead 
to the onset of a Helmholtz instability. Figure 25 represents the 
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Figure 23. Run 12, Variation of T With Respect to x/b 0 
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Figure 24. Run 12, Discrete Vortices Leading to Instabilities 
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Figure 25. Run 12, Variation of f With Respect to x/b Q 
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corresponding value of circulation for each free surface vortex at the 
time associated with Figure 24. The sharp variation in circulation 
between the vortices is noted. At the actual point of breakdown, the 
difference in circulation effectively approaches infinite values. The 
purpose of attempting to smooth the variation of circulation is to delay 
the inevitable formation of the Helmholtz instability. 

B. APPLICATION 

From Longuet-Higgins [Ref 12], a five -point smoothing routine was 
developed and implemented to smooth the value of the circulation 
when a set of adjustable parameters have been established and a set of 
specific conditions have been met. With a minor modification of the 
code, the operator can elect to either activate the smoothing subrou- 
tine continuously or require that a specified condition be met prior to 
commencing the smoothing routine. Additionally, the operator may 
specify the exact time where the routine commences looking for the 
desired conditions. In this particular investigation, the smoothing 
routine was called when the difference in circulation values between 
two adjacent vortices was greater than 0.3 and the nondimensional 
time was at least -2.1. 

Other variables that the operator can adjust include: 

1. Start time. The operator has the ability to specifically state 
when, in nondimensional time, the smoothing routine is called. 
Clearly it serves no useful purpose to smooth at early times 
because the free surface has not started deforming. Since 
smoothing is an undesirable solution of the problem in any event, 
its use should be avoided for as long as possible. A nondimen- 
sional time of -2.1 was chosen as providing a comfortable margin 
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of smoothing yet not unnecessarily introducing artificiality into 
the problem. 

2. Frequency of smoothing. Once the actual start time of smooth- 
ing has been reached, the operator may direct how often future 
smoothing attempts will be made. Arbitrarily, it was decided that 
once smoothing started, it would occur every 25 time steps. 

3. Duration of smoothing. Once smoothing has occurred for a spe- 
cific time step, the operator may (with a minor modification to 
the code) elect to automatically smooth the next “N” time steps 
prior to starting the smoothing counter once again. This selec- 
tion is useful in that the Longuet-Higgins method is not a domi- 
nating technique, and the circulation variation may be so large 
that it requires several consecutive smoothing time steps to 
remove the sharp difference which initiated the routine. In this 
application only, “N” represents the number of consecutive time 
steps that smoothing will be completed. 

4. X/b 0 position. It is unnecessary to smooth the variation of the 
circulation at points along the free surface where instabilities do 
not exist or normally form. In general, the instabilities form at an 
x/b 0 range of (0.5 < x/b 0 < 1.25) at or very close to the formation 
of the scar. 

Figure 26 represents a duplication of run 12 depicted in Figure 

23 except with the Longuet-Higgins continuous smoothing routine 

activated. For this particular run, the following smoothing routine 

selections were made: 

start time: —2.1 

smoothing duration: five time steps 

time delay 25 time steps 

x/bo range 0.5 < x/b 0 < 1.25 

Notice that the sharp peak has been reduced yet the basic form of the 

free surface deformation has been retained. 
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Figure 26. Run 12, Smoothed Variation of T With Respect to x/b Q 
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IV. RE-DISCRETIZATION 



A INTRODUCTION 

After close scrutiny of each run during this parametric study, it 
was determined that developing a method of re-discretizing the free 
surface would be useful. This task required the development of a 
routine within the model that would add a discrete vortex at a specific 
location. The critical problem would be to determine exactly when 
and where to add a free surface vortex element and what 
corresponding circulation to assign to it. The routine must have the 
ability to interrupt program execution during any specific time step, 
make the determination that a new vortex is required, and place it 
while still maintaining the continuity of the program and the general 
form of the free surface. 

B. APPLICATION 

A subroutine was developed within the code and assigned an 
activation option. If the “add vortex” option is selected, the program 
will add a discrete free surface vortex element under specific 
conditions. This routine calculates the complex distances between 
each free surface vortex element and stores the information in an 
array. When the difference of the rate of change of distance between 
two adjacent vortices exceeded an adjustable percentage value, the 
subroutine “ADDVOR” was called. For the purposes of this 
investigation, the specific value utilized was 33 percent. 
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An additional restriction was imposed upon the routine requiring 
that the two vortices under consideration be increasing the distance 
between them instead of closing the distance between them. If 
indeed the distance between two surface vortices was closing at an 
acceptable rate (i.e., greater than the 33 percent requirement), the 
subroutine would add vortices at those positions along the free surface. 
This is obviously undesirable because the vortices are already closely 
packed. Figure 27 represents a typical situation where a vortex would 
be added. 

Operator considerations include: 

1. Closing rate. As discussed earlier, the rate of change of distance 
between vortices will dictate when the program will add a vortex. 
There is no restriction as to the percentage that is chosen but it 
was determined that an acceptable range was between 30 
percent and 50 percent. 

2. Circulation. The value of the circulation of the added vortex was 
calculated algebraically between the two vortices so that the slope 
of the circulation with respect to x/bQ would be preserved 
between the (now) three surface vortices. The values of the new 
circulations are given by the following expressions: 



r; = (5*r, -r 2 )/6 



(5) 



f, = (5T 2 - r, )/6 



( 6 ) 



r„,„ =(r, + r 2 )/3 



(7) 
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Figure 27. Typical Free Surface Condition Requiring 
an Additional Vortex 
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r x represents the circulation of vortex #1 and r 2 represents the 

circulation of vortex #2 prior to the addition of the new vortex. 
/ / 

Tj and T 2 represent the new calculated vortex circulation 
values, respectively. r new represents the value of the circulation 
of the new vortex. 

3. New vortex position. In order to accurately place the new vortex, 
it was decided to calculate the new Z position utilizing a 
conservation of circulation principle. Accordingly, the position of 
the new vortex (Znew) is given by the expression: 

Znew = ((r 1 *Zj + r 2 *z 2 )-(r; *z, + r 2 *z 2 ))/r ncw (8) 



Although several placement and circulation algorithms were 
developed and tested, the ones listed in equations (5) through 
(8) performed the most adequately. Figure 28 represents the 
execution of this routine immediately following the addition of a 
vortex. This figure follows Figure 27 in time. Notice that the 
routine placed the vortex in a satisfactory manner. Figure 29 
represents the circulation strength assigned to the added vortex. 

4. Smoothing As can be seen from Figure 29, although the 
circulation is correctly assigned to the new vortex, there is a 
slight difference between the new value and the circulation values 
of its adjacent vortices. Accordingly, the operator may elect to 
immediately smooth the circulation values after adding a vortex. 
This may not be necessary, however, if the value of the circulation 
difference exceeds 0.3 of its neighbors as mentioned earlier. 
Figure 30 depicts a smoothed version of the circulation strength 
of the added vortex. 
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Figure 28. Free Surface Deformation Immediately After 

Adding a Vortex 
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Figure 29. Circulation Value of an Added Vortex 
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Figure 30. Smoothed Circulation Value of an Added Vortex 
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V. STRATIFIED FLOW EXPERIMENTS 



A EXPERIMENTAL APPARATUS 

Specific details of the experimental apparatus were discussed 
extensively by Leeker (Ref 11]. 

Experiments were conducted in a water basin approximately 
twelve feet long, three feet wide, and four feet deep. Two counter- 
rotating vortices were generated through the use of two opposing 
plates activated by a system of weights and pulleys. Flow visualization 
was performed through the use of a fluorescene dye injected through 
openings on the underside of each plate. A video camera, a VHS video 
cassette recorder, a monitor, and a date/timer were utilized as sup- 
plementary equipment. 

R PROCEDURES 

Basic experimental procedures were thoroughly discussed by 
Leeker [Ref 11]. The procedures in this investigation, although fun- 
damentally similar, differed from those of Leeker with respect to the 
addition of a sharp salt-water/fresh-water density interface. 

To obtain a sufficiently sharp density difference, 50 pounds of 
ordinary salt were thoroughly mixed with 465.7 gallons of fresh water 
in the tank bottom and allowed to settle overnight. 

From Striftos [Ref 14], the concentration is given by: 

weight [water] 

1/C= 1 + — r — : 

weight [salt] (qj 
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The density of the salt water layer is given by the expression: 

p s =p 0 *(l+.71*C) ( i o) 

Where p 0 is the density of the fresh water (62.43 lbs/ft 3 ). 

The density difference then is simply given by (Ps - Po). For the 
purposes of these experiments, C = .0127 lbs. salt/lbs. water, and 
p s = 62.993 lbs. /ft 3 . Finally, Ap =.563 lbs./ft 3 . 

A freshwater dye combination was introduced by gravity feed from 
a reservoir through small tubes which spilled onto a floating board in 
the tank. The dyed freshwater trickled over the edge of the board and 
settled on top of the salt-water layer. The tank was filled in this man- 
ner for several hours, resulting in a sharp density interface. The 
remainder of the experiments were carried out in the normal method. 
Efforts were concentrated on recording the interaction of the vortex 
pair with the interface. 

C RESULTS OF THE EXPERIMENTAL RUNS 

Approximately ten experimental runs were completed. Two runs 
(04 02 and 04 08) were chosen to discuss as representative samples. 
Analysis procedures of all experimental runs were identical to those 
discussed by Leeker [Ref 11]. 

Figure 31 depicts the right vortex (and image) center’s path 
for run 04 02. Figure 32 depicts the nondimensional velocity profile 
of the right vortex center also for run 04 02. Based on the maximum 
velocities displayed by Figure 32, the Froude number for the vortex in 
saltwater = 0.78 and for freshwater = 1.05. 
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Figure 31. Vortex Center Path for Run 04 02, F'r = 1.05 
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Figure 32. Vortex Center Nondimensional Velocity 
for Run 04 02, Fr = 1.05 
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From Figure 31, the vortex center clearly “saw" the density 
interface as a makeshift free surface and began to deflect. Upon intru- 
sion into the lower-density fluid, the velocity increased and the vortex 
path became more vertical; in effect, the vortex no longer “saw" a 
barrier. As the vortex center is influenced by the actual free surface, it 
behaved much like previous investigations [Refs 10, 11], The two dis- 
tinctly different maximum velocity peaks on Figure 32 support this 
observation. 

The vortex in run 04 08 was only visible while in the lower 
salt-water layer. Figure 33 depicts the vortex and image path. Figure 
34 represents the nondimensional velocity profile of the vortex cen- 
ter. In this run, although the center was not tracked through the 
interface, the vortex also appeared to “see” that there was a slight free 
surface above it, as evidenced by the outward slant of its path. 
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Figure 33. Vortex Center Path for Run 04 08, Fr = .38 
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Figure 34. Vortex Center Nondimensional Velocity 
for Run 04 08, Fr = .38 
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VI. CONCLUSIONS 



The numerical analysis has shown that the evolution of the free 
surface can be calculated through the use of the line vortices. No 
filtering is necessary provided that the complex velocity function is 
desingularized in a suitable manner. The analysis becomes more 
robust as the Froude number is increased. Even then, the inception of 
small-wavelength Helmholtz instability cannot be avoided. However, 
the time at which such instabilities occur can be delayed by a suitable 
arrangement of the initial position of the vortices. Fortunately, 
however, the critical time at which the instability manifests itself is 
near the time at which the wave breaking is likely to occur. This 
analysis is not equipped to deal with wave breaking. For smaller 
Froude numbers (say Fr = 0.50), the agreement between the calculated 
and measured free surface shapes is quite encouraging. A better 
agreement would have been a source of puzzlement rather than a 
source of gratification because the vortices in the experiments become 
quickly turbulent, diffuse over a large area and lose their strength due 
to the cancellation of oppositely signed vorticity in the overlapping 
regions of vortices and due to the negative vorticity generated at the 
free surface. Although the numerical model has the capability to vary 
the vortex strength with time (e.g., using a diffusing Oseen vortex), 
tills was not done. 
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For any of the Froude numbers encountered in the analysis or 
experiments, no wave train was observed on either side of the scars. 
This conclusion was somewhat anticipated on the basis of previous 
calculations. For very small Froude numbers (say, Fr < 0.15), the scars 
were very small, the vortices followed the rigid-wall path, and the 
scars were slaved to the vortices, as found previously by Sarpkaya and 
Henderson [9], 

The numerical model can be improved, within the limitations of 
the inviscid flow theory, by using vortex sheet panels rather than line 
vortices. It is also possible to use an adaptive grid system together 
with a suitable finite-difference method to solve the equations of 
motion directly. Nevertheless, the difficulties stemming from the 
occurrence of Helmholtz instability during the roll-up of the vortex 
sheets, axial instability of the vortex pair, three-dimensional 
instabilities on the free surface and the turbulent diffusion of the 
vortices remain as the most serious impediments towards an exact 
solution of the problem. 

In summary, the principal objectives of this investigation have 
been accomplished. The extensive parametric study has shown that 
the numerical model is relatively insensitive to a variety of changing 
parameters and that the Helmholtz instabilities are inevitable in the 
execution of this code beyond some critical time. 

The supplementary objective of analyzing the interaction of a 
vortex pair with a sharp density interface was also accomplished. 
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Although only a preliminary to major investigations in the future, this 
“first look” of the interaction of two counter-rotating vortices with a 
sharp density interface revealed interesting information, such as the 
considerable increase in velocity of the vortex pair after its intrusion 
into the fresh-water layer. 
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